function [Jacobian, J1, J2, J3] = cal_Jacobian_six(q, t1, t2, t3, DH_para)

Alpha = DH_para.Alpha;
A = DH_para.A;
D = DH_para.D;
Theta = DH_para.Theta;

T01 = TransformationMatrixMDH(Alpha(1), A(1), D(1), Theta(1), q(1));
T12 = TransformationMatrixMDH(Alpha(2), A(2), D(2), Theta(2), q(2));
T2_l1 = TransformationMatrixMDH(0, t1, 0, 0, 0);
T02 = T01 * T12;
T1_l1 = T12 * T2_l1;

v1 = T01(1:3, 1:3) * cross([0; 0; 1], T1_l1(1:3, 4));
v2 = T02(1:3, 1:3) * cross([0; 0; 1], T2_l1(1:3, 4));
J1 = zeros(3, 6);
J1(:, 1:2) = [v1, v2];

T01 = TransformationMatrixMDH(Alpha(1), A(1), D(1), Theta(1), q(1));
T12 = TransformationMatrixMDH(Alpha(2), A(2), D(2), Theta(2), q(2));
T23 = TransformationMatrixMDH(Alpha(3), A(3), D(3), Theta(3), q(3));
T3_l2 = TransformationMatrixMDH(Alpha(4), A(4), t2, 0, 0);
T02 = T01 * T12;
T03 = T02 * T23;
T1_l2 = T12 * T23 * T3_l2;
T2_l2 = T23 * T3_l2;
v1 = T01(1:3, 1:3) * cross([0; 0; 1], T1_l2(1:3, 4));
v2 = T02(1:3, 1:3) * cross([0; 0; 1], T2_l2(1:3, 4));
v3 = T03(1:3, 1:3) * cross([0; 0; 1], T3_l2(1:3, 4));
J2 = zeros(3, 6);
J2(:, 1:3) = [v1, v2, v3];

T01 = TransformationMatrixMDH(Alpha(1), A(1), D(1), Theta(1), q(1));
T12 = TransformationMatrixMDH(Alpha(2), A(2), D(2), Theta(2), q(2));
T23 = TransformationMatrixMDH(Alpha(3), A(3), D(3), Theta(3), q(3));
T34 = TransformationMatrixMDH(Alpha(4), A(4), D(4), Theta(4), q(4));
T45 = TransformationMatrixMDH(Alpha(5), A(5), D(5), Theta(5), q(5));
T5_l3 = TransformationMatrixMDH(Alpha(6), 0, t3, 0, 0);
T02 = T01 * T12;
T03 = T02 * T23;
T04 = T03 * T34;
T05 = T04 * T45;
T1_l3 = T12 * T23 * T34 * T45 * T5_l3;
T2_l3 = T23 * T34 * T45 * T5_l3;
T3_l3 = T34 * T45 * T5_l3;
T4_l3 = T45 * T5_l3;
v1 = T01(1:3, 1:3) * cross([0; 0; 1], T1_l3(1:3, 4));
v2 = T02(1:3, 1:3) * cross([0; 0; 1], T2_l3(1:3, 4));
v3 = T03(1:3, 1:3) * cross([0; 0; 1], T3_l3(1:3, 4));
v4 = T04(1:3, 1:3) * cross([0; 0; 1], T4_l3(1:3, 4));
v5 = T05(1:3, 1:3) * cross([0; 0; 1], T5_l3(1:3, 4));
J3 = zeros(3, 6);
J3(:, 1:5) = [v1, v2, v3, v4, v5];

T01 = TransformationMatrixMDH(Alpha(1), A(1), D(1), Theta(1), q(1));
T12 = TransformationMatrixMDH(Alpha(2), A(2), D(2), Theta(2), q(2));
T23 = TransformationMatrixMDH(Alpha(3), A(3), D(3), Theta(3), q(3));
T34 = TransformationMatrixMDH(Alpha(4), A(4), D(4), Theta(4), q(4));
T45 = TransformationMatrixMDH(Alpha(5), A(5), D(5), Theta(5), q(5));
T56 = TransformationMatrixMDH(Alpha(6), A(6), D(6), Theta(6), q(6));
T02 = T01 * T12;
T03 = T02 * T23;
T04 = T03 * T34;
T05 = T04 * T45;
T06 = T05 * T56;

T16 = T12 * T23 * T34 * T45 * T56;
T26 = T23 * T34 * T45 * T56;
T36 = T34 * T45 * T56;
T46 = T45 * T56;

v1 = T01(1:3, 1:3) * cross([0; 0; 1], T16(1:3, 4));
v2 = T02(1:3, 1:3) * cross([0; 0; 1], T26(1:3, 4));
v3 = T03(1:3, 1:3) * cross([0; 0; 1], T36(1:3, 4));
v4 = T04(1:3, 1:3) * cross([0; 0; 1], T46(1:3, 4));
v5 = T05(1:3, 1:3) * cross([0; 0; 1], T56(1:3, 4));
v6 = [0; 0; 0];

Jacobian = [v1, v2, v3, v4, v5, v6];